GitHub

Getting Started

Tensorial.jl lets you write small tensor calculations directly in Julia. This page starts with vectors and matrices, then moves on to symmetry, automatic differentiation, and block-structured states.

On this page:

Create small tensors

Vec and Mat are aliases for Tensorial vector and matrix types. If you know StaticArrays.jl, @Vec and @Mat are analogous to @SVector and @SMatrix, but they create Tensorial tensor types. Use @Tensor for higher-order tensors.

julia> using Tensorial
julia> x = @Vec [1.0, 2.0, 3.0]3-element Vec{3, Float64}:
 1.0
 2.0
 3.0
julia> A = @Mat [1.0 2.0 0.0; 0.0 3.0 4.0; 5.0 0.0 6.0]3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 1.0  2.0  0.0
 0.0  3.0  4.0
 5.0  0.0  6.0
julia> A isa AbstractArraytrue
julia> C = @Tensor rand(3, 3, 3);
julia> typeof(C)Tensor{Tuple{3, 3, 3}, Float64, 3, 27}

These tensors work with ordinary Julia array code and support tensor-specific operations:

Write tensor operations

julia> A ⊡ x3-element Vec{3, Float64}:
  5.0
 18.0
 23.0
julia> A ⊡ x ≈ A * xtrue
julia> A ⊡₂ A ≈ A ⋅ Atrue
julia> x ⊗ x3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 1.0  2.0  3.0
 2.0  4.0  6.0
 3.0  6.0  9.0

The comparisons show the familiar array equivalents. The operators themselves express tensor notation directly:

  • ⊡ contracts one index. In this example, it is the matrix-vector product.
  • ⊡₂ contracts two indices. For matrices, it is the Frobenius inner product.
  • ⊗ forms a tensor product.
Typing Unicode operators

In the Julia REPL or editors with Julia tab completion, type the LaTeX-like name and press <TAB>. ASCII aliases are also available if you prefer plain names.

OperatorTypeAlias
⊡\boxdot<TAB>contract1
⊡₂\boxdot<TAB> then \_2<TAB>contract2
⊗\otimes<TAB>tensor

Write indexed formulas with @einsum

Use @einsum when an expression is clearer in index notation. Repeated indices are summed, and indices that appear once become the free indices of the result:

julia> (@einsum A[i,k] * A[k,j]) ≈ A * Atrue
julia> (@einsum A[i,j] * A[i,j]) ≈ A ⋅ Atrue

In the first expression, k appears twice, so it is contracted. The indices i and j appear once, so they are the free indices of the result. When free indices are not written explicitly, @einsum infers them in the order they appear from left to right.

You can also give the free indices explicitly. This is useful when you want a particular output order:

julia> (@einsum (j,i) -> A[i,k] * A[k,j]) ≈ transpose(A * A)true

For a named result, put the free indices on the left-hand side:

julia> @einsum B[i,j] := A[i,k] * A[k,j];
julia> B ≈ A * Atrue

Symmetric tensors

Use symmetric(...) to compute the symmetric part of a second-order tensor and return a symmetric tensor type:

julia> ε = symmetric(@Mat [0.02 0.01 0.0; 0.01 0.00 0.0; 0.0 0.0 -0.01])3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 0.02  0.01   0.0
 0.01  0.0    0.0
 0.0   0.0   -0.01
julia> ε isa SymmetricSecondOrderTensor{3}true
julia> Tuple(ε)(0.02, 0.01, 0.0, 0.0, 0.0, -0.01)

The displayed matrix is 3×3, but the stored tuple contains only the six independent components. The symmetry is part of the tensor type, not just a property of the displayed values.

For the general @Symmetry notation and tensor-space aliases, see Tensor Types and Spaces.

Differentiate tensor formulas

Automatic differentiation is performed in the tensor space of the input. For a scalar strain-energy function, the gradient is a stress tensor and the Hessian is a tangent stiffness tensor:

julia> K = 10.0; # bulk modulus
julia> G = 5.0;  # shear modulus
julia> ψ(ε) = K/2 * tr(ε)^2 + G * (dev(ε) ⊡₂ dev(ε))ψ (generic function with 1 method)
julia> σ = gradient(ψ, ε);
julia> σ isa SymmetricSecondOrderTensor{3}true
julia> ℂ = hessian(ψ, ε);
julia> ℂ isa SymmetricFourthOrderTensor{3}true

The derivatives stay in the corresponding tensor spaces. No Voigt conversion is needed in this example.

gradient(..., :all)

Passing :all returns both the derivative and the function value. For example, g, y = gradient(f, x, :all) gives g = gradient(f, x) and y = f(x).

Direct sum (pack and unpack)

A direct sum lets you treat several blocks as one state while keeping the block structure explicit. pack builds the state, and unpack retrieves its blocks:

julia> state = pack(σ, 0.0)2-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> unpack(state, 1)3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 0.266667  0.1         0.0
 0.1       0.0666667   0.0
 0.0       0.0        -0.0333333
julia> unpack(state, 2)0.0
Infix form

⊕ is an infix spelling of pack: a ⊕ b is equivalent to pack(a, b). Type it as \oplus<TAB>.

Advanced example: block-structured return mapping

As a final, slightly more advanced example, combine the tools above in the active plastic branch of a small isotropic J2 return-mapping update. The state contains an updated stress tensor σ and a plastic multiplier increment Δγ.

The residual keeps the radial-return form explicit. Tensorial differentiates a residual with tensor and scalar blocks directly, so the Newton correction can be written in block form.

σᵗʳ = SymmetricSecondOrderTensor{3}((2.0, 0.4, 0.2, 1.2, 0.1, 0.9))
G = 5.0   # shear modulus
σy0 = 0.6 # initial yield stress
H = 2.0   # isotropic hardening modulus

q(σ) = sqrt(3/2) * norm(dev(σ)) # von Mises stress
yield_function(σ, Δγ) = q(σ) - (σy0 + H * Δγ)

function residual(x)
    σ, Δγ = unpack(x)
    # flow direction and yield-function value
    n, f = gradient(σ -> yield_function(σ, Δγ), σ, :all)
    R_σ = σ - σᵗʳ + 2G * Δγ * n
    R_γ = f
    pack(R_σ, R_γ)
end

With the residual defined, one Newton update gives a new packed state. Unpack it immediately to recover the updated variables:

x = pack(σᵗʳ, 0.0)
J = gradient(residual, x)
δx = -J \ residual(x)
xnew = x + δx
σ, Δγ = unpack(xnew)

Here x is the current Newton state, J is the Jacobian of the packed residual, δx is the Newton correction, and xnew is the updated state.

The updated stress and plastic multiplier are:

σ
3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 1.70625   0.214474   0.107237
 0.214474  1.2773     0.0536184
 0.107237  0.0536184  1.11645
Δγ
0.039112415533373614

Packing those updated blocks again gives a direct-sum state whose residual is zero for this return-mapping update:

norm(residual(σ ⊕ Δγ))
2.982917421429737e-16

For the full direct-sum explanation, including the Jacobian blocks, see Direct Sum.

Where to go next